// $Id: detector.c 479 2014-04-30 22:05:16Z alan.watson.f $

// This file is part of the control system for the Oaxaca instrument.

// This file contains the relatively high-level driver for the detectors.

////////////////////////////////////////////////////////////////////////////////

#define _POSIX_C_SOURCE 200112L

#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <time.h>

#include <sys/stat.h>
#include <unistd.h>

#include "oaxaca-config.h"
#include "oaxaca.h"
#include "ao.h"
#include "array.h"
#include "detector.h"
#include "ids.h"
#include "mprintf.h"
#include "rfits.h"
#include "spots.h"
#include "wfits.h"
#include "wpng.h"

////////////////////////////////////////////////////////////////////////////////

#include "minmax.h"

////////////////////////////////////////////////////////////////////////////////

static array_short_t *mean = NULL;
static array_short_t *dark = NULL;

////////////////////////////////////////////////////////////////////////////////

double 
oaxaca_detector_scale(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->scale;
  else
    return 0.0;
}

void 
oaxaca_detector_set_scale(int index, double scale)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->scale = scale;
}

double 
oaxaca_detector_pixelfrequency(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->pixelfrequency;
  else
    return 0.0;
}

void 
oaxaca_detector_set_pixelfrequency(int index, double pixelfrequency)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->pixelfrequency = pixelfrequency;
}

int 
oaxaca_detector_window_cx(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->cx;
  else
    return 0;
}

void 
oaxaca_detector_set_window_cx(int index, int cx)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->cx = cx;
}

int 
oaxaca_detector_window_cy(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->cy;
  else
    return 0;
}

void 
oaxaca_detector_set_window_cy(int index, int cy)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->cy = cy;
}

int 
oaxaca_detector_window_sx(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->window.sx;
  else
    return 0;
}

void 
oaxaca_detector_set_window_sx(int index, int sx)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->window.sx = sx;
}

int 
oaxaca_detector_window_sy(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->window.sy;
  else
    return 0;
}

void 
oaxaca_detector_set_window_sy(int index, int sy)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->window.sy = sy;
}

int 
oaxaca_detector_window_nx(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->window.nx;
  else
    return 0;
}

void 
oaxaca_detector_set_window_nx(int index, int nx)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->window.nx = nx;
}

int 
oaxaca_detector_window_ny(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->window.ny;
  else
    return 0;
}

void 
oaxaca_detector_set_window_ny(int index, int ny)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->window.ny = ny;
}

double 
oaxaca_detector_minframefrequency(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->minframefrequency;
  else
    return 0.0;
}

double 
oaxaca_detector_maxframefrequency(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->maxframefrequency;
  else
    return 0.0;
}

double 
oaxaca_detector_framefrequency(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->framefrequency;
  else
    return 0.0;
}

void 
oaxaca_detector_set_framefrequency(int index, double framefrequency)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->framefrequency = framefrequency;
}

double 
oaxaca_detector_actualframefrequency(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->actualframefrequency;
  else
    return 0.0;
}

double 
oaxaca_detector_frametime(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->frametime;
  else
    return 0.0;
}

double 
oaxaca_detector_exposuretime(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->exposuretime;
  else
    return 0.0;
}

void 
oaxaca_detector_set_exposuretime(int index, double exposuretime)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->exposuretime = exposuretime;
}

double 
oaxaca_detector_integrationtime(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->integrationtime;
  else
    return 0.0;
}

void 
oaxaca_detector_set_integrationtime(int index, double integrationtime)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    detector->integrationtime = integrationtime;
}

unsigned long
oaxaca_detector_framesperintegration(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->framesperintegration;
  else
    return 0;
}

unsigned long
oaxaca_detector_integrationsperexposure(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->integrationsperexposure;
  else
    return 0;
}

unsigned long
oaxaca_detector_integration(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->integration;
  else
    return 0;
}

unsigned long
oaxaca_detector_frame(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->frame % detector->framesperintegration;
  else
    return 0;
}

double
oaxaca_detector_zmin(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->zmin;
  else
    return 0;
}

double
oaxaca_detector_zmax(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->zmax;
  else
    return 0;
}

double
oaxaca_detector_zmean(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->zmean;
  else
    return 0;
}

double
oaxaca_detector_fwhm(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->fwhm;
  else
    return 0;
}

double
oaxaca_detector_eshm(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->eshm;
  else
    return 0;
}

double
oaxaca_detector_spatialsigma(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->spatialsigma;
  else
    return 0;
}

double
oaxaca_detector_temporalsigma(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (oaxaca_shm != NULL && detector->initialized)
    return detector->temporalsigma;
  else
    return 0;
}

void
oaxaca_detector_set_abort(int index, int abort)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  detector->abort = abort;
}

unsigned long
oaxaca_detector_search_radius(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->searchradius;
}

void
oaxaca_detector_set_search_radius(int index, unsigned long radius)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  detector->searchradius = radius;
}

unsigned long
oaxaca_detector_refinement_radius(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->refinementradius;
}

void
oaxaca_detector_set_refinement_radius(int index, unsigned long radius)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  detector->refinementradius = radius;
}

unsigned long
oaxaca_detector_ns(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->ns;
}

void
oaxaca_detector_set_ns(int index, unsigned long ns)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  detector->ns = ns;
}

double
oaxaca_detector_s_vector_ref(int index, unsigned long i)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->s[i];
}

double
oaxaca_detector_r_vector_ref(int index, unsigned long i)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->r[i];
}

void
oaxaca_detector_r_vector_set(int index, unsigned long i, double r)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  detector->r[i] = r;
}

double
oaxaca_detector_ds_mean_vector_ref(int index, unsigned long i)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->dsmean[i];
}

double
oaxaca_detector_ds_sigma_vector_ref(int index, unsigned long i)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->dssigma[i];
}

double
oaxaca_detector_i_matrix_ref(int index, unsigned long i, unsigned long j)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (i < OAXACA_NA_MAX && j < OAXACA_NS_MAX)
    return detector->i_matrix[i][j];
  else
    return 0.0;
}

void
oaxaca_detector_i_matrix_set(int index, unsigned long i, unsigned long j, double x)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (i < OAXACA_NA_MAX && j < OAXACA_NS_MAX)
    detector->i_matrix[i][j] = x;
  else
    oaxaca_error("bad i-matrix indices.");
}

double
oaxaca_detector_c_matrix_ref(int index, unsigned long i, unsigned long j)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (i < OAXACA_NA_MAX && j < OAXACA_NS_MAX)
    return detector->c_matrix[i][j];
  else
    return 0.0;
}

void
oaxaca_detector_c_matrix_set(int index, unsigned long i, unsigned long j, double x)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  if (i < OAXACA_NA_MAX && j < OAXACA_NS_MAX)
    detector->c_matrix[i][j] = x;
  else
    oaxaca_error("bad c-matrix indices.");
}

double
oaxaca_detector_g_vector_ref(int index, unsigned long i)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  return detector->g[i];
}

void
oaxaca_detector_g_vector_set(int index, unsigned long i, double g)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];
  detector->g[i] = g;
}

////////////////////////////////////////////////////////////////////////////////

void
oaxaca_detector_open(int index, int id, int ns)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  oaxaca_info("opening detector.");
  detector->initialized = 0;
  oaxaca_shm->detector[index] = oaxaca_ids_initialize(id);
  detector->scale = 1.0;
  detector->ns = ns;
  for (int i = 0; i < ns; ++i) {
    detector->r[i] = 0.0;
    detector->s[i] = 0.0;
  }
  for (unsigned long i = 0; i < OAXACA_NA_MAX; ++i)
    for (unsigned long j = 0; j < OAXACA_NS_MAX; ++j)
      detector->i_matrix[i][j] = 0.0;
  for (unsigned long i = 0; i < OAXACA_NA_MAX; ++i)
    for (unsigned long j = 0; j < OAXACA_NS_MAX; ++j)
      detector->c_matrix[i][j] = 0.0;
  for (unsigned long i = 0; i < OAXACA_NA_MAX; ++i)
    detector->g[i] = 0.0;
  detector->initialized = 1;
  oaxaca_info("finished opening detector.");
}

void
oaxaca_detector_configure(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  detector->zmin                 = 0.0;
  detector->zmax                 = 0.0;
  detector->zmean                = 0.0;
  detector->fwhm                 = 0.0;
  detector->eshm                 = 0.0;
  detector->spatialsigma         = 0.0;
  detector->temporalsigma        = 0.0;
  detector->actualframefrequency = 0.0;
  detector->totaldelay           = 0.0;
  oaxaca_ao_clear(index);
  
  oaxaca_info("configuring detector.");
  oaxaca_shm->detector[index] = oaxaca_ids_configure(oaxaca_shm->detector[index]);    
  oaxaca_info("finished configuring detector.");

  oaxaca_info("reading the dark image.");
  free_array_short(dark);
  dark = NULL;
  FILE *fp = fopen("dark.fits", "rb");
  if (fp != NULL) {
    dark = oaxaca_rfits(fp);
    fclose(fp);
  }
  if (dark == NULL) {
    oaxaca_info("unable to read the dark image.");
  } else if (dark->nx != detector->window.nx || 
      dark->ny != detector->window.ny) {
    oaxaca_info("the dark image does not match the size of the window.");
    free_array_short(dark);
    dark = NULL;
  } else {
    oaxaca_info("finished reading the dark image.");
  }
}

void
oaxaca_detector_close(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  oaxaca_info("closing detector.");
  if (oaxaca_shm != NULL)
    detector->initialized = 0;
  oaxaca_info("finished closing detector.");
}

////////////////////////////////////////////////////////////////////////////////

static void
integrations_log(const char *filename, int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  FILE *fp = fopen("integrations.csv", "a");
  if (fp == NULL)
    oaxaca_error("unable to open integrations log file.");

  fprintf(fp, "\"%s\"", filename);
  
  time_t seconds = time(NULL);
  struct tm tm = *gmtime(&seconds);
  char yyyymmdd[] = "YYYY-MM-DD";
  strftime(yyyymmdd, sizeof(yyyymmdd), "%Y-%m-%d", &tm);
  char hhmmss[] = "hh:mm:ss";
  strftime(hhmmss, sizeof(hhmmss), "%H:%M:%S", &tm);
  fprintf(fp, ",\"%s %s\"", yyyymmdd, hhmmss);
  fprintf(fp, ",\"%s\"", yyyymmdd);
  fprintf(fp, ",\"%s\"", hhmmss);

  fprintf(fp, ",%.3f", detector->fwhm);
  fprintf(fp, ",%.3f", detector->fwhm * detector->scale);
  fprintf(fp, ",%.3f", detector->hartmannfwhm);
  fprintf(fp, ",%.3f", detector->hartmannfwhm * detector->scale);
  fprintf(fp, ",%.3f", detector->scale);
  fprintf(fp, ",%.3e", detector->eshm);
  fprintf(fp, ",%.3e", detector->framefrequency);
  fprintf(fp, ",%.3e", detector->actualframefrequency);
  fprintf(fp, ",%.6f", detector->frametime);
  fprintf(fp, ",%.6f", detector->integrationtime);
  fprintf(fp, ",%lu", detector->integration);
  fprintf(fp, ",%.3e", detector->zmin);
  fprintf(fp, ",%.3e", detector->zmax);
  fprintf(fp, ",%.3e", detector->zmean);
  fprintf(fp, ",%.3e", detector->zmean);
  fprintf(fp, ",%.3e", detector->spatialsigma);
  fprintf(fp, ",%.3e", detector->temporalsigma);
  
  fprintf(fp, ",\"%c\"", oaxaca_shm->wfs.command);

  volatile oaxaca_detector_t *wfs_detector = &oaxaca_shm->detector[OAXACA_DETECTOR_INDEX_WFS];  
  fprintf(fp, ",%.3f", wfs_detector->fwhm);
  fprintf(fp, ",%.3f", wfs_detector->fwhm * detector->scale);
  fprintf(fp, ",%.3f", wfs_detector->hartmannfwhm);
  fprintf(fp, ",%.3f", wfs_detector->hartmannfwhm * detector->scale);
  fprintf(fp, ",%.3e", wfs_detector->eshm);
  fprintf(fp, ",%.3f", wfs_detector->scale);
  fprintf(fp, ",%.3e", wfs_detector->framefrequency);
  fprintf(fp, ",%.3e", wfs_detector->actualframefrequency);
  fprintf(fp, ",%.6f", wfs_detector->frametime);
  fprintf(fp, ",%.6f", wfs_detector->integrationtime);
  fprintf(fp, ",%lu", wfs_detector->integration);
  fprintf(fp, ",%.3e", wfs_detector->zmin);
  fprintf(fp, ",%.3e", wfs_detector->zmax);
  fprintf(fp, ",%.3e", wfs_detector->zmean);
  fprintf(fp, ",%.3e", wfs_detector->zmean);
  fprintf(fp, ",%.3e", wfs_detector->spatialsigma);
  fprintf(fp, ",%.3e", wfs_detector->temporalsigma);
  for (int i = 0; i < OAXACA_NA_MAX; ++i)  
    fprintf(fp, ",%.3f", wfs_detector->g[i]);

  for (int i = 0; i < OAXACA_NA_MAX; ++i)  
    fprintf(fp, ",%.3f", oaxaca_shm->wfc.amean[i]);
  for (int i = 0; i < OAXACA_NA_MAX; ++i)  
    fprintf(fp, ",%.3f", oaxaca_shm->wfc.asigma[i]);
  for (int i = 0; i < OAXACA_NA_MAX; ++i)  
    fprintf(fp, ",%.3f", oaxaca_shm->wfc.damean[i]);
  for (int i = 0; i < OAXACA_NA_MAX; ++i)  
    fprintf(fp, ",%.3f", oaxaca_shm->wfc.dasigma[i]);
  for (int i = 0; i < OAXACA_NA_MAX; ++i)  
    fprintf(fp, ",%.3f", oaxaca_shm->wfc.asaturated[i]);

  fprintf(fp, ",%.5f", oaxaca_shm->telescope.alpha);
  fprintf(fp, ",%.5f", oaxaca_shm->telescope.delta);
  fprintf(fp, ",%.5f", oaxaca_shm->telescope.ha);
  fprintf(fp, ",%.5f", oaxaca_shm->telescope.zd);

  fprintf(fp, ",%.3e", oaxaca_shm->dimm.lambda);
  fprintf(fp, ",%.3f", oaxaca_shm->dimm.r0);
  fprintf(fp, ",%.3f", oaxaca_shm->dimm.fwhm);
  fprintf(fp, ",%.3f", oaxaca_shm->dimm.fwhm / (detector->fwhm * detector->scale));

  fprintf(fp, "\n");
  
  fclose(fp);  
}

////////////////////////////////////////////////////////////////////////////////

static const char *
get_basename(void)
{
  struct timeval tv;
  gettimeofday(&tv, NULL);
  time_t seconds = tv.tv_sec;
  int mseconds = tv.tv_usec / 1000;
  struct tm tm = *gmtime(&seconds);
  char yyyymmdd[] = "YYYYMMDD";
  strftime(yyyymmdd, sizeof(yyyymmdd), "%Y%m%d", &tm);
  char hhmmss[] = "hhmmss";
  strftime(hhmmss, sizeof(hhmmss), "%H%M%S", &tm);
  return mprintf("%sT%s.%03d", yyyymmdd, hhmmss, mseconds);
}

static void
write_fits_file(const char *filename, int index, array_short_t *a)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  FILE *fp = fopen(filename, "wb");
  if (fp == NULL)
    oaxaca_error("unable to open FITS file.");

  oaxaca_wfits_write_basic_header(fp, a->nx, a->ny);
  oaxaca_wfits_write_date(fp, "DATE");
  oaxaca_wfits_write_string(fp, "ORIGNAME", filename);
  oaxaca_wfits_write_double(fp, "PXFR", detector->pixelfrequency, "%.2e");
  oaxaca_wfits_write_long(fp, "WNSX", detector->window.sx);
  oaxaca_wfits_write_long(fp, "WNSY", detector->window.sy);
  oaxaca_wfits_write_long(fp, "WNNX", detector->window.nx);
  oaxaca_wfits_write_long(fp, "WNNY", detector->window.ny);
  oaxaca_wfits_write_double(fp, "FRFR", detector->framefrequency, "%.2f");
  oaxaca_wfits_write_double(fp, "ACFRFR", detector->actualframefrequency, "%.2f");
  oaxaca_wfits_write_double(fp, "FRTM", detector->frametime, "%.6f");
  oaxaca_wfits_write_double(fp, "INTM", detector->integrationtime, "%.6f");
  oaxaca_wfits_write_double(fp, "FRPERIN", detector->framesperintegration, "%.0f");

  oaxaca_wfits_write_data(fp, a);

  if (fclose(fp) != 0)
    oaxaca_error("while writing FITS file.");
}

////////////////////////////////////////////////////////////////////////////////

static void (*frame_call_back) (int, array_short_t *) = NULL;
static void (*integration_call_back) (int, array_short_t *) = NULL;

static array_double_t *s1 = NULL;
static array_double_t *s2 = NULL;

static void
integration_start(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  oaxaca_info("starting integration %ld.", detector->integration);
  
  size_t nx = detector->window.nx;
  size_t ny = detector->window.ny;

  if (s1 == NULL || s1->nx != nx || s1->ny != ny) {
    free_array_double(s1);
    s1 = alloc_array_double(ny, nx);
  }
  for (int iy = 0; iy < ny; ++iy)
    for (int ix = 0; ix < nx; ++ix)
      s1->z[iy][ix] = 0.0;    

  if (s2 == NULL || s2->nx != nx || s2->ny != ny) {
    free_array_double(s2);
    s2 = alloc_array_double(ny, nx);  
  }
  for (int iy = 0; iy < ny; ++iy)
    for (int ix = 0; ix < nx; ++ix)
      s2->z[iy][ix] = 0.0;
}

static void
integration_handler(int index, array_short_t *frame)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  size_t nx = detector->window.nx;
  size_t ny = detector->window.ny;

  for (int iy = 0; iy < ny; ++iy)
    for (int ix = 0; ix < nx; ++ix) {
      double z = (double) frame->z[iy][ix];
      s1->z[iy][ix] += z;
      s2->z[iy][ix] += z * z;
    }
}

static void
integration_end(int index)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  long integration = detector->integration;
  
  oaxaca_info("finished integration %ld.", integration);
  
  size_t nx = detector->window.nx;
  size_t ny = detector->window.ny;
  
  double s0 = detector->framesperintegration;

  // Determine the minimum, maximum, the mean, the total signal, the spatial noise, the FWHM, and the ESHM.
  {
    double zmin = +INFINITY;
    double zmax = -INFINITY;
    double t0 = nx * ny;
    double t1 = 0.0;
    double t2 = 0.0;
    for (int iy = 0; iy < ny; ++iy)
      for (int ix = 0; ix < nx; ++ix) {
        double z = s1->z[iy][ix] / s0;
        zmin = min(z, zmin);
        zmax = max(z, zmax);
        t1 += z;
        t2 += z * z;
      }
    double mean = t1 / t0;
    double sigma = sqrt((t0 * t2 - t1 * t1) / t0 / (t0 - 1));
    detector->zmin = zmin;
    detector->zmax = zmax;
    detector->zmean = mean;
    detector->spatialsigma = sigma;
  }
  
  // Determine the mean total signal for each spot.
  int ns = detector->ns;
  double zsum = 0.0;
  for (int i = 0; i < ns / 2; ++i) {
    int rx = (int) (detector->r[2 * i + 0] + 0.5);
    int ry = (int) (detector->r[2 * i + 1] + 0.5);
    int searchradius = detector->searchradius;
    int sx = imax(rx - searchradius, 0);
    int sy = imax(ry - searchradius, 0);
    int ex = imin(rx + searchradius, nx - 1);
    int ey = imin(ry + searchradius, ny - 1);
    for (int iy = sy; iy <= ey; ++iy) {
      for (int ix = sx; ix <= ex; ++ix) {
        if ((ix - rx) * (ix - rx) + (iy - ry) * (iy - ry) <= searchradius * searchradius) {
          double z = s1->z[iy][ix] / s0;
          zsum += z;
        }
      }
    }
  }
  double meantotalsignal = zsum / (ns / 2);
  oaxaca_debug("total signal is %.0f.", meantotalsignal);

  // Determine the FWHM and ESHM for each spot.
  double sumsqfwhm = 0.0;
  double sumeshm = 0.0;
  for (int i = 0; i < ns / 2; ++i) {
    int rx = (int) (detector->r[2 * i + 0] + 0.5);
    int ry = (int) (detector->r[2 * i + 1] + 0.5);
    int searchradius = detector->searchradius;
    int sx = imax(rx - searchradius, 0);
    int sy = imax(ry - searchradius, 0);
    int ex = imin(rx + searchradius, nx - 1);
    int ey = imin(ry + searchradius, ny - 1);
    double zmax = -INFINITY;
    for (int iy = sy; iy <= ey; ++iy) {
      for (int ix = sx; ix <= ex; ++ix) {
        if ((ix - rx) * (ix - rx) + (iy - ry) * (iy - ry) <= searchradius * searchradius) {
          double z = s1->z[iy][ix] / s0;
          zmax = max(z, zmax);
        }
      }
    }
    double nhm  = 0.0;
    double eshm = 0.0;
    for (int iy = sy; iy <= ey; ++iy) {
      for (int ix = sx; ix <= ex; ++ix) {
        if ((ix - rx) * (ix - rx) + (iy - ry) * (iy - ry) <= searchradius * searchradius) {
          double z = s1->z[iy][ix] / s0;
          if (z > 0.5 * zmax) {
            nhm  += 1.0;
            eshm += z;
          }
        }
      }
    }
    double pi = 4.0 * atan(1.0);
    double fwhm = 2.0 * sqrt(nhm / pi);
    oaxaca_debug("FWHM of spot %d is %.1f pixels.", i, fwhm);
    sumsqfwhm += fwhm * fwhm;
    sumeshm += eshm;
  }
  if (ns > 0) {
    double fwhm = sqrt(sumsqfwhm / (ns / 2));
    double eshm = sumeshm / (ns / 2);
    oaxaca_info("spot FWHM is %.2f arcsec.", fwhm * detector->scale);
    oaxaca_debug("spot ESHM is %.0f.", eshm);
    detector->fwhm = fwhm;
    detector->eshm = eshm;
  }
  
  // Determine the mean temporal noise.
  if (s0 == 1) {
    detector->temporalsigma = 0.0;
  } else {
    double t0 = nx * ny;
    double t2 = 0.0;
    for (int iy = 0; iy < ny; ++iy)
      for (int ix = 0; ix < nx; ++ix) {
        double r = sqrt((s0 * s2->z[iy][ix] - s1->z[iy][ix] * s1->z[iy][ix]) / s0 / (s0 - 1));
        t2 += r * r;
      }
    double sigma = sqrt(t2 / t0);
    detector->temporalsigma = sigma / sqrt(s0);
  }

  // Determine the mean.
  if (mean == NULL || mean->nx != nx || mean->ny != ny) {
    free_array_short(mean);
    mean = alloc_array_short(ny, nx);
  }
  for (int iy = 0; iy < ny; ++iy)
    for (int ix = 0; ix < nx; ++ix)
      mean->z[iy][ix] = (short) round(s1->z[iy][ix] / s0);

  const char *basename = get_basename();
  char *fits_filename = mprintf("%s.fits", basename);
  integrations_log(fits_filename, index);

  signal(SIGCHLD, SIG_IGN);
  if (fork() == 0) {
  
    oaxaca_info("writing FITS file \"%s\".", fits_filename);
    write_fits_file(fits_filename, index, mean);

    char *tmp_filename = mprintf("%s.tmp", basename);
    unlink(tmp_filename);
    if (link(fits_filename, tmp_filename) != 0 ||
        rename(tmp_filename, "latest.fits") != 0)
      oaxaca_error("unable to create latest.fits.");
    
    int ns = detector->ns;
    double s[ns];
    for (int i = 0; i < ns; ++i)
      s[i] = detector->r[i];
    unsigned long searchradius = detector->searchradius;
    unsigned long refinementradius = detector->refinementradius;
    oaxaca_wpng_write_file(mean, ns, s, searchradius, refinementradius);

    free(fits_filename);
    free(tmp_filename);
    free((void *) basename);
    
    _exit(0);
  }
  
  if (integration_call_back != NULL)
    integration_call_back(index, mean);

}

static void
correct_bias_and_background(array_short_t *frame)
{
  // Frames have bias instability which manifests itself as horizontal
  // striping modulated by the column parity. We attempt to remove this
  // by subtracting the mean in the odd and even columns in the first
  // 10% and last 10% of each row from the corresponding columns of each
  // row. This also removes the background. We then repeat for the columns.
  
  int nx = frame->nx;
  int ny = frame->ny;

  for (int iy = 0; iy < ny; ++iy) {
    long n[2] = { 0, 0 };
    long sum[2] = { 0, 0 };
    for (int ix = 0; ix < nx / 10; ++ix) {
      int i = ix & 1;
      ++n[i];
      sum[i] += frame->z[iy][ix];
    }
    for (int ix = 9 * nx / 10; ix < nx; ++ix) {
      int i = ix & 1;
      ++n[i];
      sum[i] += frame->z[iy][ix];
    }
    short mean[2] = { sum[0] / n[0], sum[1] / n[1] };
    for (int ix = 0; ix < nx; ++ix) {
      int i = ix & 1;
      frame->z[iy][ix] -= mean[i];
    }
  }

  for (int ix = 0; ix < nx; ++ix) {
    long n = 0;
    long sum = 0;
    for (int iy = 0; iy < ny / 10; ++iy) {
      ++n;
      sum += frame->z[iy][ix];
    }
    for (int iy = 9 * ny / 10; iy < ny; ++iy) {
      ++n;
      sum += frame->z[iy][ix];
    }
    short mean = sum / n;
    for (int iy = 0; iy < ny; ++iy) {
      frame->z[iy][ix] -= mean;
    }
  }

}

static int
frame_handler(int index, double actualframefrequency, array_short_t *frame)
{
  volatile oaxaca_detector_t *detector = &oaxaca_shm->detector[index];

  if (detector->abort) {
    oaxaca_info("aborting.");
    detector->abort = 0;
    return 0;
  }
  
  struct timeval tv;
  gettimeofday(&tv, NULL);
  double timestamp = tv.tv_sec + tv.tv_usec * 1e-6;
  detector->timestamp = timestamp;
  
  detector->actualframefrequency = actualframefrequency;

  if (detector->frame % detector->framesperintegration == 0)
    integration_start(index);

  size_t nx = detector->window.nx;
  size_t ny = detector->window.ny;

  for (int iy = 0; iy < ny; ++iy)
    for (int ix = 0; ix < nx; ++ix)
      frame->z[iy][ix] <<= 2;
  if (dark != NULL) {
    for (int iy = 0; iy < ny; ++iy)
      for (int ix = 0; ix < nx; ++ix)
        frame->z[iy][ix] -= dark->z[iy][ix];
    correct_bias_and_background(frame);
  }

  if (frame_call_back != NULL)
    frame_call_back(index, frame);

  integration_handler(index, frame);

  ++detector->frame;

  if (detector->frame % detector->framesperintegration == 0) {
    integration_end(index);
    ++detector->integration;
    if (detector->integration == detector->integrationsperexposure)
      return 0;
  }

  return 1;
}

static void
take_frames(int index)
{
  oaxaca_detector_t detector = oaxaca_shm->detector[index];

  detector.zmin          = 0.0;
  detector.zmax          = 0.0;
  detector.zmean         = 0.0;
  detector.fwhm          = 0.0;
  detector.eshm          = 0.0;
  detector.spatialsigma  = 0.0;
  detector.temporalsigma = 0.0;
  detector.rmsdsmean     = 0.0;
  detector.rmsdssigma    = 0.0;
  detector.hartmannfwhm  = 0.0;

  detector.actualframefrequency = 0.0;

  oaxaca_debug("frame time is %.6fs.", detector.frametime);

  oaxaca_debug("requested integration time is %.6fs.", detector.integrationtime);
  if (detector.integrationtime == 0)
    detector.framesperintegration = 0;
  else {
    detector.framesperintegration = round(detector.integrationtime / detector.frametime);
    if (detector.framesperintegration == 0)
      detector.framesperintegration = 1;
  }
  detector.integrationtime = detector.framesperintegration * detector.frametime;
  detector.frame = 0;
  oaxaca_debug("actual integration time is %.6fs.", detector.integrationtime);
  oaxaca_debug("frames per integration is %lu.", detector.framesperintegration);
  
  oaxaca_debug("requested exposure time is %.6fs.", detector.exposuretime);
  if (detector.exposuretime == 0)
    detector.integrationsperexposure = 0;
  else {
    detector.integrationsperexposure = round(detector.exposuretime / detector.integrationtime);
    if (detector.integrationsperexposure == 0)
      detector.integrationsperexposure = 1;
  }
  detector.exposuretime = detector.integrationsperexposure * detector.integrationtime;
  detector.integration = 0;
  oaxaca_debug("actual exposure time is %.6fs.", detector.exposuretime);
  oaxaca_debug("integrations per exposure is %lu.", detector.integrationsperexposure);
  
  oaxaca_shm->detector[index] = detector;

  oaxaca_ids_run(index, frame_handler);
}

////////////////////////////////////////////////////////////////////////////////

void
oaxaca_detector_expose(int index)
{
  oaxaca_info("starting exposure.");
  oaxaca_ao_clear(index);
  take_frames(index);
  oaxaca_info("finished exposure.");
}


void
oaxaca_detector_ao(int index)
{
  oaxaca_info("starting correcting.");
  oaxaca_ao_clear(index);
  frame_call_back       = oaxaca_ao_frame_call_back;
  integration_call_back = oaxaca_ao_integration_call_back;
  take_frames(index);
  frame_call_back       = NULL;
  integration_call_back = NULL;
  oaxaca_info("finished correcting.");
}

////////////////////////////////////////////////////////////////////////////////

void
oaxaca_detector_dark(int index)
{
  oaxaca_info("starting dark exposure.");
  oaxaca_ao_clear(index);
  free_array_short(dark);
  dark = NULL;
  take_frames(index);
  dark = mean;
  mean = NULL;
  oaxaca_info("finished dark exposure.");
  write_fits_file("dark.fits", index, dark);
}

////////////////////////////////////////////////////////////////////////////////